/*
** FloatSymmetricLevinson.cs
**
** Class that implements UDL factorisation of the inverse of a symmetric square
** Toeplitz matrix using the Levinson Algorithm.
**
** Copyright (c) 2003-2004, dnAnalytics Project. All rights reserved.
*/

#region Using directives

using System;
using dnA.Exceptions;

#endregion Using directives

namespace dnA.Math
{

	/// <summary>
	/// A Levinson solver for symmetric square Toeplitz systems of <c>float</c> type.
	/// </summary>
	/// <threadsafety static="true" instance="true"/>
	/// <remarks>
	/// This class provides members for inverting a symmetric square Toeplitz matrix
	/// (see <see cref="GetInverse"/> member), calculating the determinant of the matrix
	/// (see <see cref="GetDeterminant"/> member) and solving linear systems associated
	/// with the matrix (see <see cref="Solve"/> members).
	/// <para>
	/// The class implements an <B>UDL</B> decomposition of the inverse of the Toeplitz matrix.
	/// The decomposition is based upon Levinson's algorithm. As a consequence, all operations
	/// require approximately <B>N</B> squared FLOPS, where <B>N</B> is the matrix order. This
	/// is significantly faster than Cholesky's factorization for symmetric matrices, which
	/// requires <B>N</B> cubed FLOPS.
	/// </para>
	/// <para>
	/// A requirement of Levinson's algorithm is that all the leading sub-matrices and the principal
	/// matrix must be non-singular. During the decomposition, sub-matrices and the principal matrix are
	/// checked to ensure that they are non-singular. When a singular matrix is found, the decomposition
	/// is halted and an internal flag is set. The <see cref="IsSingular"/> property may be used to access
	/// the flag, to determine if any singular matrices were detected.
	/// </para>
	/// <para>
	/// It has been shown that Levinson's algorithm is weakly numerically stable for positive-definite
	/// Toeplitz matrices. It is usual to restrict the use of the algorithm to such
	/// matrix types. The <see cref="IsPositiveDefinite"/> property may be checked to verify that the
	/// Toeplitz matrix is positive-definite.
	/// </para>
	/// <para>
	/// If one of the leading sub-matrices or the principal matrix is near-singular, then the accuracy
	/// of the decomposition will be degraded. An estimate of the resulting error is provided and may be
	/// accessed with the <see cref="Error"/> property. This estimate is only valid for positive-definite
	/// matrices.
	/// </para>
	/// <para>
	/// A outline of this approach to the UDL decomposition of inverse Toeplitz matrices is
	/// found in the following reference:
	/// </para>
	/// <para>
	/// <EM>Sun-Yuan Kung and Yu Hen Hu</EM>, A Highly Concurrent Algorithm and Pipelined
	/// Architecture for Solving Toeplitz Systems, IEEE Transactions on Acoustics,
	/// Speech and Signal Processing, Volume ASSP-31, Number 1, Febuary 1983, pages
	/// 66 - 75.
	/// </para>
	/// </remarks>
	/// <example>
	/// The following simple example illustrates the use of the class:
	/// <para>
	/// <code escaped="true">
	/// using System;
	/// using dnA.Exceptions;
	/// using dnA.Math;
	/// using System.IO;
	///
	/// namespace Example_3
	/// {
	///
	///   class Application
	///   {
	///
	///     // format string for matrix/vector elements
	///     private const string frmStr = " 0.000E+000;-0.000E+000";
	///
	///     // application entry point
	///     [STAThread]
	///     public static void Main(string[] args)
	///     {
	///
	///       // create Levinson solver
	///       FloatSymmetricLevinson fsl = new FloatSymmetricLevinson(1.0f, 0.5f, 0.2f);
	///
	///       // display the Toeplitz matrix
	///       FloatMatrix T = fsl.GetMatrix();
	///       Console.WriteLine("Matrix:: {0} ", T.ToString(frmStr));
	///       Console.WriteLine();
	///
	///       // check if matrix is singular
	///       Console.WriteLine("Singular:               {0}", fsl.IsSingular);
	///
	///       // check if matrix is positive definite
	///       Console.WriteLine("Positive Definite:      {0}", fsl.IsPositiveDefinite);
	///
	///       // get error for inverse
	///       Console.WriteLine("Cybenko Error Estimate: {0}", fsl.Error.ToString("E3"));
	///
	///       // get the determinant
	///       Console.WriteLine("Determinant:            {0}", fsl.GetDeterminant().ToString("E3"));
	///       Console.WriteLine();
	///
	///       // get the inverse
	///       FloatMatrix Inv = fsl.GetInverse();
	///       Console.WriteLine("Inverse:: {0} ", Inv.ToString(frmStr));
	///       Console.WriteLine();
	///
	///       // solve a linear system
	///       FloatVector X = fsl.Solve(4.0f, -1.0f, 3.0f);
	///       Console.WriteLine("X:: {0} ", X.ToString(frmStr));
	///       Console.WriteLine();
	///
	///     }
	///
	///   }
	///
	/// }
	/// </code>
	/// </para>
	/// <para>
	/// The application generates the following results:
	/// </para>
	/// <para>
	/// <code escaped="true">
	/// Matrix:: rows: 3, cols: 3
	///  1.000E+000,  5.000E-001,  2.000E-001
	///  5.000E-001,  1.000E+000,  5.000E-001
	///  2.000E-001,  5.000E-001,  1.000E+000
	///
	/// Singular:               False
	/// Positive Definite:      True
	/// Cybenko Error Estimate: 4.087E-007
	/// Determinant:            5.600E-001
	///
	/// Inverse:: rows: 3, cols: 3
	///  1.339E+000, -7.143E-001,  8.929E-002
	/// -7.143E-001,  1.714E+000, -7.143E-001
	///  8.929E-002, -7.143E-001,  1.339E+000
	///
	/// X:: Length: 3
	///  6.339E+000, -6.714E+000,  5.089E+000
	/// </code>
	/// </para>
	/// </example>
	sealed public class FloatSymmetricLevinson : Algorithm
	{

		#region Constants

		/// <summary>
		/// The unit roundoff for float type.
		/// </summary>
		/// <remarks>
		/// It is assumed that the <see cref="float"/> type has a 23-bit mantissa.
		/// </remarks>
		private const float UNITROUNDOFF = 1.1920929e-7F;

		#endregion Constants

		#region Fields

		/// <summary>
		/// The left-most column of the Toeplitz matrix.
		/// </summary>
		private readonly FloatVector m_LeftColumn;

		/// <summary>
		/// The system order.
		/// </summary>
		private readonly int m_Order;

		/// <summary>
		/// lower triangular matrix.
		/// </summary>
		private float[][] m_LowerTriangle;

		/// <summary>
		/// The diagonal vector.
		/// </summary>
		private float[] m_Diagonal;

		/// <summary>
		/// Flag to indicate if system is singular.
		/// </summary>
		private bool m_IsSingular;

		/// <summary>
		/// Flag to indicate if system is positive definite.
		/// </summary>
		private bool m_IsPositiveDefinite;

		/// <summary>
		/// The Cybenko error-bound
		/// </summary>
		private float m_CybenkoError;

		#endregion Fields

		#region Properties

		/// <summary>
		/// Get the order of the Toeplitz matrix.
		/// </summary>
		public int Order
		{
			get
			{
				return m_Order;
			}
		}

		/// <summary>
		/// Check if the Toeplitz matrix or any leading sub-matrices are singular.
		/// </summary>
		/// <remarks>
		/// If the Toeplitz matrix or any leading sub-matrices are singular, it is
		/// not possible to complete the <b>UDL</b> decomposition using the Levinson
		/// algorithm.
		/// </remarks>
		public bool IsSingular
		{
			get
			{
				// make sure the factorisation is completed
				Compute();
				// return the flag value
				return m_IsSingular;
			}
		}

		/// <summary>
		/// Check if the Toeplitz matrix is positive definite.
		/// </summary>
		/// <remarks>
		/// It has only been shown that the Levinson algorithm is weakly numerically stable
		/// for symmetric positive-definite Toeplitz matrices.
		/// Based on empirical results, it appears that the Levinson algorithm
		/// gives reasonable accuracy for many symmetric indefinite matrices.
		/// It may be desirable to restrict the use of this class to positive-definite matrices
		/// to guarantee accuracy.
		/// </remarks>
		public bool IsPositiveDefinite
		{
			get
			{
				// make sure the factorisation is completed
				Compute();
				// return the flag value
				return m_IsPositiveDefinite;
			}
		}

		///<summary>
		/// Get the lower triangle matrix of the UDL factorisation.
		/// </summary>
		/// <remarks>
		/// It is recommended that the <see cref="IsSingular"/> property
		/// be checked to see if the decomposition was completed, before attempting
		/// to obtain the lower triangle matrix.
		/// </remarks>
		/// <exception cref="SingularMatrixException">
		/// The Toeplitz matrix or one of the the leading sub-matrices is singular.
		/// </exception>
		public FloatMatrix L
		{
			get
			{
				// make sure the factorisation is completed
				Compute();

				// check if there was a result
				if (m_IsSingular == true)
				{
					throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
				}

				// copy jagged array into a FloatMatrix
				FloatMatrix Lower = new FloatMatrix(m_Order);
				for (int i = 0; i < m_Order; i++ )
				{
#if MANAGED
					// managed implementation
					m_LowerTriangle[i].CopyTo(Lower.data[i], 0);
#else
					// native implementation
					for( int j = 0; j < i+1; j++ )
					{
						Lower[i, j] = m_LowerTriangle[i][j];
					}
#endif
				}

				return Lower;
			}
		}

		///<summary>
		/// Get the diagonal matrix of the UDL factorisation.</summary>
		/// <remarks>
		/// It is recommended that the <see cref="IsSingular"/> property
		/// be checked to see if the decomposition was completed, before attempting
		/// to obtain the diagonal matrix.
		/// </remarks>
		/// <exception cref="SingularMatrixException">
		/// The Toeplitz matrix or one of the the leading sub-matrices is singular.
		/// </exception>
		public FloatMatrix D
		{
			get
			{
				// make sure the factorisation is completed
				Compute();

				// check if there was a result
				if (m_IsSingular == true)
				{
					throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
				}

				// copy diagonal vector into a FloatMatrix
				FloatMatrix Diagonal = new FloatMatrix(m_Order);
				for (int i = 0; i < m_Order; i++ )
				{
#if MANAGED
					// managed implementation
					Diagonal.data[i][i] = m_Diagonal[i];
#else
					// native implementation
					Diagonal.data[i*m_Order+i] = m_Diagonal[i];
#endif
				}

				return Diagonal;
			}
		}

		///<summary>
		/// Get the upper triangle matrix of the UDL factorisation.</summary>
		/// <remarks>
		/// It is recommended that the <see cref="IsSingular"/> property
		/// be checked to see if the decomposition was completed, before attempting
		/// to obtain the upper triangle matrix.
		/// </remarks>
		/// <exception cref="SingularMatrixException">
		/// The Toeplitz matrix or one of the the leading sub-matrices is singular.
		/// </exception>
		public FloatMatrix U
		{
			get
			{
				// make sure the factorisation is completed
				Compute();

				// check if there was a result
				if (m_IsSingular == true)
				{
					throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
				}

				// copy jagged array into a FloatMatrix and then calculate the transpose
				FloatMatrix Upper = this.L.GetTranspose();

				return Upper;
			}
		}

		/// <summary>
		/// Get an error estimate for the inverse matrix.
		/// </summary>
		/// <remarks>
		/// This estimate is approximate and is only valid for positive definite matrices.
		/// It is useful for checking the accurary of the <b>UDL</b> decomposition. If the
		/// Toeplitz matrix or one of the leading sub-matrices is near-singular,
		/// the accuracy of the decomposition will be degraded. This error estimate
		/// can be used to identify such situations.
		/// <para>
		/// The estimate is based on the Cybenko error-bound for the inverse
		/// and is a relative error. When calculating this error it is assumed that 
		/// the <see cref="float"/> type utilises a 23-bit mantissa.
		/// </para>
		/// <para>
		/// The error estimate will range in value between 1.19e-7 and <see cref="float.PositiveInfinity"/>.
		/// </para>
		/// </remarks>
		/// <example>
		/// The following example calculates the inverse of a symmetric positive-definite Toeplitz matrix.
		/// An exact inverse was obtained previously using a computer algebra system. The exact inverse is
		/// compared to the inverse calculated by this class. The actual error is identified and compared
		/// to the Cybenko error-bound.
		/// <para>
		/// <code>
		/// using System;
		/// using dnA.Exceptions;
		/// using dnA.Math;
		/// using System.IO;
		///
		/// namespace Example_4
		/// {
		///
		///   class Application
		///   {
		///
		///     // The main entry point for the application.
		///     [STAThread]
		///     public static void Main(string[] args)
		///     {
		///       // exact inverse from yacas
		///       FloatMatrix exact = new FloatMatrix(5);
		///       exact[4, 4] = exact[0, 0] = +1.3577842E+000f;
		///       exact[3, 4] = exact[4, 3] = exact[1, 0] = exact[0, 1] = -5.9039646E-001f;
		///       exact[2, 4] = exact[4, 2] = exact[2, 0] = exact[0, 2] = -1.0830325E-001f;
		///       exact[1, 4] = exact[4, 1] = exact[3, 0] = exact[0, 3] = -5.9423022E-002f;
		///       exact[4, 0] = exact[0, 4] = -5.8145106E-002f;
		///       exact[3, 3] = exact[1, 1] = +1.6120124E+000f;
		///       exact[3, 2] = exact[2, 3] = exact[2, 1] = exact[1, 2] = -5.4584837E-001f;
		///       exact[3, 1] = exact[1, 3] = -8.7102652E-002f;
		///       exact[2, 2] = +1.6180506E+000f;
		///
		///       // create levinson solver
		///       FloatSymmetricLevinson fsl = new FloatSymmetricLevinson(1.0f, 1.0f/2.0f, 1.0f/3.0f, 1.0f/4.0f, 1.0f/5.0f);
		///
		///       // identify relative error
		///       FloatMatrix Error = fsl.GetInverse() - exact;
		///       float e = Error.GetL1Norm() / exact.GetL1Norm();
		///
		///       // display results
		///       Console.WriteLine("Observed Error:         {0}", e.ToString("E3"));
		///       Console.WriteLine("Cybenko Error Estimate: {0}", fsl.Error.ToString("E3"));
		///       Console.WriteLine();
		///     }
		///
		///   }
		///
		/// }
		/// </code>
		/// </para>
		/// <para>
		/// The application generates the following results.
		/// </para>
		/// <para>
		/// <code escaped="true">
		/// Observed Error:         6.110E-008
		/// Cybenko Error Estimate: 5.520E-007
		/// </code>
		/// </para>
		/// </example>
		public float Error
		{
			get
			{
				// make sure the factorisation is completed
				Compute();
				// return the error estimate
				return m_CybenkoError;
			}
		}

		#endregion Properties

		#region Constructors

		/// <overloads>
		/// There are two permuations of the constructor, both require a parameter corresponding
		/// to the left-most column of a Toeplitz matrix.
		/// </overloads>
		/// <summary>
		/// Constructor with <c>FloatVector</c> parameter.
		/// </summary>
		/// <param name="T">The left-most column of the Toeplitz matrix.</param>
		/// <exception cref="ArgumentNullException">
		/// <B>T</B> is a null reference.
		/// </exception>
		/// <exception cref="RankException">
		/// The length of <B>T</B> is zero.
		/// </exception>
		public FloatSymmetricLevinson(FloatVector T)
		{
			// check parameter
			if (T == null)
			{
				throw new System.ArgumentNullException("T");
			}
			else if (T.Length == 0)
			{
				throw new RankException("The length of T is zero.");
			}

			// save the vector
			m_LeftColumn = T.Clone();
			m_Order = m_LeftColumn.Length;

			// allocate memory for lower triangular matrix
			m_LowerTriangle = new float[m_Order][];
			for (int i = 0; i < m_Order; i++)
			{
				m_LowerTriangle[i] = new float[i+1];
			}

			// allocate memory for diagonal
			m_Diagonal = new float[m_Order];
		}

		/// <summary>
		/// Constructor with <c>float</c> array parameter.
		/// </summary>
		/// <param name="T">The left-most column of the Toeplitz matrix.</param>
		/// <exception cref="ArgumentNullException">
		/// <B>T</B> is a null reference.
		/// </exception>
		/// <exception cref="RankException">
		/// The length of <B>T</B> is zero.
		/// </exception>
		public FloatSymmetricLevinson(params float[] T)
		{
			// check parameter
			if (T == null)
			{
				throw new System.ArgumentNullException("T");
			}
			else if (T.Length == 0)
			{
				throw new RankException("The length of T is zero.");
			}

			// save the vector
			m_LeftColumn = new FloatVector(T);
			m_Order = m_LeftColumn.Length;

			// allocate memory for lower triangular matrix
			m_LowerTriangle = new float[m_Order][];
			for (int i = 0; i < m_Order; i++)
			{
				m_LowerTriangle[i] = new float[i+1];
			}

			// allocate memory for diagonal
			m_Diagonal = new float[m_Order];
		}

		#endregion Constructors

		#region Protected Members

		/// <summary>
		/// Calculate the UDL decomposition of the inverse Toeplitz matrix.
		/// </summary>
		/// <remarks>
		/// The member checks each iteration for a singular principal sub-matrix.
		/// If a singular sub-matrix is detected, the private boolean <EM>m_IsSingular</EM>
		/// is set to <B>true</B>, the <EM>m_IsPositiveDefinite</EM>EM> flag is set <B>false</B> and the
		/// decomposition halted.
		/// <para>
		/// When the decomposition is completed, the lower triangular matrix is stored
		/// in the jagged array <EM>m_LowerTriangle</EM>, the diagonal elements
		/// in the array <EM>m_Diagonal</EM> and the private boolean <EM>m_IsSingular</EM>
		/// is set to <B>false</B>.
		/// </para>
		/// <para>
		/// The diagonal elements are checked to determine if the Toeplitz matrix is
		/// positive definite. The <EM>m_IsPositiveDefinite flag</EM>EM> is set <B>true</B>, if all
		/// diagonal elements are positive, otherwise it is set to <B>false</B>.
		/// </para>
		/// <para>
		/// During the calculation the reflection coefficients are used to calculate the
		/// Cybenko error estimate which is stored in the member <EM>m_CybenkoError</EM>.
		/// </para>
		/// </remarks>
		protected override void InternalCompute()
		{
			int i, j, l;			// index/loop variables
			float Inner;			// inner product
			float K;				// reflection coefficient
			float[] B;				// reference to previous order coefficients
			float[] A;				// reference to current order coefficients


			// check if principal diagonal is zero
			if (m_LeftColumn[0] == 0.0f)
			{
				m_IsSingular = true;
				m_IsPositiveDefinite = false;
				m_CybenkoError = float.MaxValue;
				return;
			}

			// setup zero order solution
			B = m_LowerTriangle[0];
			B[0] = 1.0f;
			m_Diagonal[0] = 1.0f / m_LeftColumn[0];
			m_CybenkoError = UNITROUNDOFF;

			// solve systems of increasing order
			for (i = 1; i < m_Order; i++, B = A)
			{
				// calculate inner product
				Inner = 0.0f;
				for ( j = 0, l = 1; j < i; j++, l++ )
				{
					Inner += m_LeftColumn[l] * B[j];
				}

				// calculate the reflection coefficient
				K = -Inner * m_Diagonal[i-1];

				// get the current low triangle row
				A = m_LowerTriangle[i];

				// update low triangle elements
				A[0] = 0.0f;
				Array.Copy(B, 0, A, 1, i);
				for (j = 0, l = i - 1; j < i; j++, l--)
				{
					A[j] += K * B[l];
				}

				// check if singular sub-matrix
				if (K*K == 1.0f)
				{
					m_IsSingular = true;
					m_IsPositiveDefinite = false;
					m_CybenkoError = float.MaxValue;
					return;
				}

				// update diagonal
				m_Diagonal[i] = m_Diagonal[i - 1]/(1.0f - K * K);

				// update error estimate
				m_CybenkoError *= (1.0f + System.Math.Abs(K)) / (1.0f - System.Math.Abs(K));
			}

			// check if system is positive-definite
			for (i = 0, m_IsPositiveDefinite = true; (i < m_Order) && (m_IsPositiveDefinite == true); i++)
			{
				if (m_Diagonal[i] < 0.0f)
				{
					m_IsPositiveDefinite = false;
				}
			}

			// solution completed
			m_IsSingular = false;
			return;
		}

		#endregion Protected Members

		#region Public Methods

		/// <summary>
		/// Get a vector that represents the left-most column of the Toplitz matrix.
		/// </summary>
		public FloatVector GetVector()
		{
			return new FloatVector(m_LeftColumn);
		}

		/// <summary>
		/// Get a copy of the Toeplitz matrix.
		/// </summary>
		public FloatMatrix GetMatrix()
		{
			int i, j;

			// allocate memory for the matrix
			FloatMatrix tm = new FloatMatrix(m_Order);

#if MANAGED
			// fill top row
			float[] top = tm.data[0];
			Array.Copy(m_LeftColumn.data, 0, top, 0, m_Order);

			if (m_Order > 1)
			{
				// fill bottom row (reverse order)
				float[] bottom = tm.data[m_Order - 1];

				for (i = 0, j = m_Order - 1; i < m_Order; i++, j--)
				{
					bottom[i] = m_LeftColumn[j];
				}

				// fill rows in-between
				for (i = 1, j = m_Order - 1 ; j > 1; i++)
				{
					Array.Copy(top, 0, tm.data[i], i, j--);
					Array.Copy(bottom, j, tm.data[i], 0, i);
				}
			}
#else
			if (m_Order > 1)
			{
				float[] top = new float[m_Order];
				Array.Copy(m_LeftColumn.data, 0, top, 0, m_Order);
				tm.SetRow(0, top);

				// fill bottom row (reverse order)
				float[] bottom = new float[m_Order];

				for (i = 0, j = m_Order - 1; i < m_Order; i++, j--)
				{
					bottom[i] = m_LeftColumn[j];
				}

				// fill rows in-between
				for (i = 1, j = m_Order - 1 ; j > 0; i++)
				{
					float[] temp = new float[m_Order];
					Array.Copy(top, 0, temp, i, j--);
					Array.Copy(bottom, j, temp, 0, i);
					tm.SetRow(i, temp);
				}
			}
			else
			{
				Array.Copy(m_LeftColumn.data, 0, tm.data, 0, m_Order);
			}
#endif

			return tm;
		}

		/// <summary>
		/// Get the determinant of the Toeplitz matrix.
		/// </summary>
		/// <returns>The determinant.</returns>
		/// <remarks>
		/// It is recommended that the <see cref="IsSingular"/> property
		/// be checked to see if the decomposition was completed, before attempting
		/// to obtain the determinant.
		/// </remarks>
		/// <exception cref="SingularMatrixException">
		/// The Toeplitz matrix or one of the the leading sub-matrices is singular.
		/// </exception>
		public float GetDeterminant()
		{
			float Determinant;

			// make sure the factorisation is completed
			Compute();

			// if any of the matrices is singular give up
			if (m_IsSingular == true)
			{
				throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
			}
			else
			{
				Determinant = 1.0f;
				for (int i = 0; i < m_Order; i++)
				{
					Determinant *= m_Diagonal[i];
				}
				Determinant = 1.0f / Determinant;
			}

			return Determinant;
		}

		/// <overloads>
		/// Solve a symmetric square Toeplitz system.
		/// </overloads>
		/// <summary>
		/// Solve a symmetric square Toeplitz system with a right-side vector.
		/// </summary>
		/// <param name="Y">The right-hand side of the system.</param>
		/// <returns>The solution vector.</returns>
		/// <exception cref="ArgumentNullException">
		/// Parameter <B>Y</B> is a null reference.
		/// </exception>
		/// <exception cref="RankException">
		/// The length of <B>Y</B> is not equal to the number of rows in the Toeplitz matrix.
		/// </exception>
		/// <exception cref="SingularMatrixException">
		/// The Toeplitz matrix or one of the the leading sub-matrices is singular.
		/// </exception>
		/// <remarks>
		/// This member solves the linear system <B>TX</B> = <B>Y</B>, where <B>T</B> is
		/// the symmetric square Toeplitz matrix, <B>X</B> is the unknown solution vector
		/// and <B>Y</B> is a known vector.
		/// <para>
		/// The class implicitly decomposes the inverse Toeplitz matrix into a <b>UDL</b> factorisation
		/// using the Levinson algorithm, and then calculates the solution vector.
		/// </para>
		/// </remarks>
		public FloatVector Solve(FloatVector Y)
		{
			FloatVector X;

			// check parameters
			if (Y == null)
			{
				throw new System.ArgumentNullException("Y");
			}
			else if (m_Order != Y.Length)
			{
				throw new RankException("The length of Y is not equal to the number of rows in the Toeplitz matrix.");
			}

			Compute();

			if (m_IsSingular == true)
			{
				throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
			}

			int i, j, l;			// index/loop variables
			float Inner;			// inner product
			float G;				// scaling constant
			float[] A;				// reference to current order coefficients

			// allocate memory for solution
			X = new FloatVector(m_Order);

			// setup zero order solution
			X[0] = Y[0] / m_LeftColumn[0];

			// solve systems of increasing order
			for (i = 1; i < m_Order; i++)
			{
				// calculate inner product
				Inner = Y[i];
				for (j = 0, l = i; j < i; j++, l--)
				{
					Inner -= X[j] * m_LeftColumn[l];
				}
				// get the current predictor coefficients row
				A = m_LowerTriangle[i];

				// update the solution vector
				G = Inner * m_Diagonal[i];
				for (j = 0; j <= i; j++)
				{
					X[j] += G * A[j];
				}

			}

			return X;
		}

		/// <overloads>
		/// Solve a symmetric square Toeplitz system.
		/// </overloads>
		/// <summary>
		/// Solve a symmetric square Toeplitz system with a right-side array.
		/// </summary>
		/// <param name="Y">The right-hand side of the system.</param>
		/// <returns>The solution vector.</returns>
		/// <exception cref="ArgumentNullException">
		/// Parameter <B>Y</B> is a null reference.
		/// </exception>
		/// <exception cref="RankException">
		/// The length of <B>Y</B> is not equal to the number of rows in the Toeplitz matrix.
		/// </exception>
		/// <exception cref="SingularMatrixException">
		/// The Toeplitz matrix or one of the the leading sub-matrices is singular.
		/// </exception>
		/// <remarks>
		/// This member solves the linear system <B>TX</B> = <B>Y</B>, where <B>T</B> is
		/// the symmetric square Toeplitz matrix, <B>X</B> is the unknown solution vector
		/// and <B>Y</B> is a known vector.
		/// <para>
		/// The class implicitly decomposes the inverse Toeplitz matrix into a <b>UDL</b> factorisation
		/// using the Levinson algorithm, and then calculates the solution vector.
		/// </para>
		/// </remarks>
		public FloatVector Solve(params float[] Y)
		{

			FloatVector X;

			// check parameters
			if (Y == null)
			{
				throw new System.ArgumentNullException("Y");
			}
			else if (m_Order != Y.Length)
			{
				throw new RankException("The length of Y is not equal to the number of rows in the Toeplitz matrix.");
			}

			Compute();

			if (m_IsSingular == true)
			{
				throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
			}

			int i, j, l;			// index/loop variables
			float Inner;			// inner product
			float G;				// scaling constant
			float[] A;				// reference to current order coefficients

			// allocate memory for solution
			X = new FloatVector(m_Order);

			// setup zero order solution
			X[0] = Y[0] / m_LeftColumn[0];

			// solve systems of increasing order
			for (i = 1; i < m_Order; i++)
			{
				// calculate inner product
				Inner = Y[i];
				for (j = 0, l = i; j < i; j++, l--)
				{
					Inner -= X[j] * m_LeftColumn[l];
				}
				// get the current predictor coefficients row
				A = m_LowerTriangle[i];

				// update the solution vector
				G = Inner * m_Diagonal[i];
				for (j = 0; j <= i; j++)
				{
					X[j] += G * A[j];
				}

			}

			return X;
		}

		/// <summary>
		/// Solve a symmetric square Toeplitz system with a right-side matrix.
		/// </summary>
		/// <param name="Y">The right-hand side of the system.</param>
		/// <returns>The solution matrix.</returns>
		/// <exception cref="ArgumentNullException">
		/// Parameter <B>Y</B> is a null reference.
		/// </exception>
		/// <exception cref="RankException">
		/// The number of rows in <B>Y</B> is not equal to the number of rows in the Toeplitz matrix.
		/// </exception>
		/// <exception cref="SingularMatrixException">
		/// The Toeplitz matrix or one of the the leading sub-matrices is singular.
		/// </exception>
		/// <remarks>
		/// This member solves the linear system <B>TX</B> = <B>Y</B>, where <B>T</B> is
		/// a symmetric square Toeplitz matrix, <B>X</B> is the unknown solution matrix
		/// and <B>Y</B> is a known matrix.
		/// <para>
		/// The class implicitly decomposes the inverse Toeplitz matrix into a <b>UDL</b> factorisation
		/// using the Levinson algorithm, and then calculates the solution matrix.
		/// </para>
		/// </remarks>
		public FloatMatrix Solve(FloatMatrix Y)
		{
			FloatMatrix X;

			// check parameters
			if (Y == null)
			{
				throw new System.ArgumentNullException("Y");
			}
			else if (m_Order != Y.ColumnLength)
			{
				throw new RankException("The numer of rows in Y is not equal to the number of rows in the Toeplitz matrix.");
			}

			Compute();

			if (m_IsSingular == true)
			{
				throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
			}

			int M = Y.RowLength;
			int i, j, l, m;			// index/loop variables
			float[] Inner;			// inner product
			float[] G;				// scaling constant
			float[] A;				// reference to current order coefficients
			float scalar;

			// allocate memory for solution
			X = new FloatMatrix(m_Order, M);
			Inner = new float[M];
			G = new float[M];

			// setup zero order solution
			scalar = 1.0f / m_LeftColumn[0];
			for (m = 0; m < M; m++)
			{
#if MANAGED
				X.data[0][m] = scalar * Y.data[0][m];
#else

				X.data[m*m_Order] = scalar * Y.data[m*m_Order];
#endif
			}

			// solve systems of increasing order
			for (i = 1; i < m_Order; i++)
			{
				// calculate inner product
				for (m = 0; m < M; m++)
				{
#if MANAGED
					Inner[m] = Y.data[i][m];
#else
					Inner[m] = Y.data[m*m_Order+i];
#endif
				}

				for (j = 0, l = i; j < i; j++, l--)
				{
					scalar = m_LeftColumn[l];
					for (m = 0; m < M; m++)
					{
#if MANAGED
						Inner[m] -= scalar * X.data[j][m];
#else
						Inner[m] -= scalar * X.data[m*m_Order+j];
#endif
					}
				}

				// get the current predictor coefficients row
				A = m_LowerTriangle[i];

				// update the solution matrix
				for (m = 0; m < M; m++)
				{
					G[m] = m_Diagonal[i] * Inner[m];
				}
				for (j = 0; j <= i; j++)
				{
					scalar = A[j];
					for (m = 0; m < M; m++)
					{
#if MANAGED
						X.data[j][m] += scalar * G[m];
#else
						X.data[m*m_Order+j] += scalar * G[m];
#endif
					}
				}
			}

			return X;
		}

		/// <summary>
		/// Get the inverse of the Toeplitz matrix.
		/// </summary>
		/// <returns>The inverse matrix.</returns>
		/// <exception cref="SingularMatrixException">
		/// The Toeplitz matrix or one of the the leading sub-matrices is singular.
		/// </exception>
		/// <remarks>
		/// The class implicitly decomposes the inverse Toeplitz matrix into a <b>UDL</b> factorisation
		/// using the Levinson algorithm, before using Trench's algorithm to complete
		/// the calculation of the inverse.
		/// <para>
		/// Trench's algorithm requires approximately <b>N</b> squared FLOPS, compared to <b>N</b> cubed FLOPS
		/// if we simply multiplied the <b>UDL</b> factors (<b>N</b> is the matrix order).
		/// </para>
		/// </remarks>
		public FloatMatrix GetInverse()
		{
			Compute();

			if (m_IsSingular == true)
			{
				throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
			}

			FloatMatrix I = new FloatMatrix(m_Order);						// the solution matrix
			float[] A = m_LowerTriangle[m_Order-1];
			float A1, A2, scale;

#if MANAGED

			float[] current, previous;										// references to rows in the solution
			int i, j, k, l;

			// setup the first row in wedge
			scale = m_Diagonal[m_Order-1];
			current = I.data[0];
			for (i = 0, j = m_Order - 1; i < m_Order; i++, j--)
			{
				current[i] = scale* A[j];
			}

			// calculate values in the rest of the wedge
			for (i = 1; i < (1 + m_Order) / 2; i++)
			{
				previous = current;
				current = I.data[i];
				A1 = A[m_Order - i - 1];
				A2 = A[i - 1];
				for (j = i, k = i - 1, l = m_Order - i - 1; j < m_Order - i; j++, k++, l--)
				{
					current[j] = previous[k] + scale * (A1 * A[l] - A2 * A[k]);
				}
			}

#else

			int i, j, k, l;

			// setup the first row in wedge
			scale = m_Diagonal[m_Order-1];
			for (i = 0, j = m_Order - 1; i < m_Order; i++, j--)
			{
				I[0, i] = scale* A[j];
			}

			// calculate values in the rest of the wedge
			for (i = 1; i < (1 + m_Order) / 2; i++)
			{
				A1 = A[m_Order - i - 1];
				A2 = A[i - 1];
				for (j = i, k = i - 1, l = m_Order - i - 1; j < m_Order - i; j++, k++, l--)
				{
					I[i, j] = I[i - 1, k] + scale * (A1 * A[l] - A2 * A[k]);
				}
			}

#endif

			// this is symmetric matrix ...
			for (i = 0; i < (1 + m_Order) / 2; i++)
			{
				for (j = i; j < m_Order - i; j++)
				{
					I[j, i] = I[i, j];
				}
			}

			// and a persymmetric matrix.
			for (i = 0, j = m_Order - 1; i < m_Order; i++, j--)
			{
				for (k = 0, l = m_Order - 1; k < j; k++, l--)
				{
					I[l, j] = I[i, k];
				}
			}

			return I;
		}

		#endregion Public Methods

		#region Public Static Methods

		/// <overloads>
		/// Solve a symmetric square Toeplitz system.
		/// </overloads>
		/// <summary>
		/// Solve a symmetric square Toeplitz system with a right-side vector.
		/// </summary>
		/// <param name="T">The left-most column of the Toeplitz matrix.</param>
		/// <param name="Y">The right-side vector of the system.</param>
		/// <returns>The solution vector.</returns>
		/// <exception cref="ArgumentNullException">
		/// <B>T</B> and/or <B>Y</B> are null references
		/// </exception>
		/// <exception cref="RankException">
		/// The length of <B>T</B> does not match the length of <B>Y</B>.
		/// </exception>
		/// <exception cref="SingularMatrixException">
		/// The Toeplitz matrix or one of the the leading sub-matrices is singular.
		/// </exception>
		/// <remarks>
		/// This method solves the linear system <B>AX</B> = <B>Y</B>. Where
		/// <B>T</B> is a symmetric square Toeplitz matrix, <B>X</B> is an unknown
		/// vector and <B>Y</B> is a known vector.
		/// <para>
		/// This static member combines the <b>UDL</b> decomposition and the calculation of the solution into a
		/// single algorithm. When compared to the non-static member it requires minimal data storage
		/// and suffers from no speed penalty.
		/// </para>
		/// </remarks>
		public static FloatVector Solve(FloatVector T, FloatVector Y)
		{

			FloatVector X;

			// check parameters
			if (T == null)
			{
				throw new System.ArgumentNullException("T");
			}
			else if (Y == null)
			{
				throw new System.ArgumentNullException("Y");
			}
			else if (T.Length != Y.Length)
			{
				throw new RankException("The length of T and Y are not equal.");
			}
			else
			{

				// allocate memory
				int N = T.Length;
				X = new FloatVector(N);                    // solution vector
				float e;                                   // prediction error

				// setup zero order solution
				e = T[0];
				if (e == 0.0f)
				{
					throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
				}
				X[0] = Y[0] / T[0];

				if (N > 1)
				{
					FloatVector a = new FloatVector(N - 1);   // prediction coefficients
					FloatVector Z = new FloatVector(N - 1);   // temporary storage vector
					float g;                                   // reflection coefficient
					float inner;                               // inner product
					float k;
					int i, j, l;

					// calculate solution for successive orders
					for (i = 1; i < N; i++)
					{

						// calculate first inner product
						inner = T[i];
						for (j = 0, l = i - 1; j < i - 1; j++, l--)
						{
							inner += a[j] * T[l];
						}

						// update predictor coefficients
						g = -(inner / e);
						for (j = 0, l = i - 2; j < i - 1; j++, l--)
						{
							Z[j] = a[j] + g * a[l];
						}

						// copy vector
						for (j = 0; j < i - 1; j++)
						{
							a[j] = Z[j];
						}

						a[i - 1] = g;

						e *= (1.0f - g * g);
						if (e == 0.0f)
						{
							throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
						}

						// calculate second inner product
						inner = Y[i];
						for (j = 0, l = i; j < i; j++, l--)
						{
							inner -= X[j] * T[l];
						}

						// update solution vector
						k = inner / e;
						for (j = 0, l = i - 1; j < i; j++, l--)
						{
							X[j] = X[j] + k * a[l];
						}
						X[j] = k;

					}

				}

			}

			return X;

		}

		/// <summary>
		/// Solve a symmetric square Toeplitz system with a right-side matrix.
		/// </summary>
		/// <param name="T">The left-most column of the Toeplitz matrix.</param>
		/// <param name="Y">The right-side matrix of the system.</param>
		/// <returns>The solution matrix.</returns>
		/// <exception cref="ArgumentNullException">
		/// <B>T</B> and/or <B>Y</B> are null references
		/// </exception>
		/// <exception cref="RankException">
		/// The length of <B>T</B> does not match the number of rows in <B>Y</B>.
		/// </exception>
		/// <exception cref="SingularMatrixException">
		/// The Toeplitz matrix or one of the the leading sub-matrices is singular.
		/// </exception>
		/// <remarks>
		/// This method solves the linear system <B>AX</B> = <B>Y</B>. Where
		/// <B>T</B> is a symmetric square Toeplitz matrix, <B>X</B> is an unknown
		/// matrix and <B>Y</B> is a known matrix.
		/// <para>
		/// This static member combines the <b>UDL</b> decomposition and the calculation of the solution into a
		/// single algorithm. When compared to the non-static member it requires minimal data storage
		/// and suffers from no speed penalty.
		/// </para>
		/// </remarks>
		public static FloatMatrix Solve(FloatVector T, FloatMatrix Y)
		{

			FloatMatrix X;

			// check parameters
			if (T == null)
			{
				throw new System.ArgumentNullException("T");
			}
			else if (Y == null)
			{
				throw new System.ArgumentNullException("Y");
			}
			else if (T.Length != Y.ColumnLength)
			{
				throw new RankException("The length of T and Y are not equal.");
			}
			else
			{

				// allocate memory
				int N = T.Length;
				int M = Y.RowLength;
				X = new FloatMatrix(N, M);                 // solution matrix
				FloatVector Z = new FloatVector(N);       // temporary storage vector
				float e;                                   // prediction error
				int i, j, l, m;

				// setup zero order solution
				e = T[0];
				if (e == 0.0f)
				{
					throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
				}
				for (m = 0; m < M; m++)
				{
					X[0, m] = Y[0,m] / T[0];
				}

				if (N > 1)
				{

					FloatVector a = new FloatVector(N - 1);   // prediction coefficients
					float p;                                   // reflection coefficient
					float inner;                               // inner product
					float k;

					// calculate solution for successive orders
					for (i = 1; i < N; i++)
					{

						// calculate first inner product
						inner = T[i];
						for (j = 0, l = i - 1; j < i - 1; j++, l--)
						{
							inner += a[j] * T[l];
						}

						// update predictor coefficients
						p = -(inner / e);
						for (j = 0, l = i - 2; j < i - 1; j++, l--)
						{
							Z[j] = a[j] + p * a[l];
						}

						// copy vector
						for (j = 0; j < i - 1; j++)
						{
							a[j] = Z[j];
						}

						a[i - 1] = p;
						e *= (1.0f - p * p);

						if (e == 0.0f)
						{
							throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
						}

						// update the solution matrix
						for (m = 0; m < M; m++)
						{

							// retrieve a copy of solution column
							for (j = 0; j < i; j++)
							{
								Z[j] = X[j, m];
							}

							// calculate second inner product
							inner = Y[i, m];
							for (j = 0, l = i; j < i; j++, l--)
							{
								inner -= Z[j] * T[l];
							}

							// update solution vector
							k = inner / e;
							for (j = 0, l = i - 1; j < i; j++, l--)
							{
								Z[j] = Z[j] + k * a[l];
							}
							Z[j] = k;

							// store solution column in matrix
							for (j = 0; j <= i; j++)
							{
								X[j, m] = Z[j];
							}

						}

					}

				}

			}

			return X;
		}

		/// <summary>
		/// Solve the Yule-Walker equations for a symmetric square Toeplitz system
		/// </summary>
		/// <param name="R">The left-most column of the Toeplitz matrix.</param>
		/// <returns>The solution vector.</returns>
		/// <exception cref="ArgumentNullException">
		/// <B>R</B> is a null reference.
		/// </exception>
		/// <exception cref="ArgumentOutOfRangeException">
		/// The length of <B>R</B> must be greater than one.
		/// </exception>
		/// <exception cref="SingularMatrixException">
		/// The Toeplitz matrix or one of the the leading sub-matrices is singular.
		/// </exception>
		/// <remarks>
		/// This member is used to solve the Yule-Walker system <B>AX</B> = -<B>a</B>,
		/// where <B>A</B> is a symmetric square Toeplitz matrix, constructed
		/// from the elements <B>R[0]</B>, ..., <B>R[N-2]</B> and
		/// the vector <B>a</B> is constructed from the elements
		/// <B>R[1]</B>, ..., <B>R[N-1]</B>.
		/// <para>
		/// Durbin's algorithm is used to solve the linear system. It requires
		/// approximately the <b>N</b> squared FLOPS to calculate the
		/// solution (<b>N</b> is the matrix order).
		/// </para>
		/// </remarks>
		public static FloatVector YuleWalker(FloatVector R)
		{

			FloatVector a;

			// check parameters
			if (R == null)
			{
				throw new System.ArgumentNullException("R");
			}
			else if (R.Length < 2)
			{
				throw new System.ArgumentOutOfRangeException("R", "The length of R must be greater than 1.");
			}
			else
			{

				int N = R.Length - 1;
				a = new FloatVector(N);                    // prediction coefficients
				FloatVector Z = new FloatVector(N);		// temporary storage vector
				float e;				                    // predictor error
				float inner;                               // inner product
				float g;                                   // reflection coefficient
				int i, j, l;

				// setup first order solution
				e = R[0];
				if (e == 0.0f)
				{
					throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
				}
				g = -R[1] / R[0];
				a[0] = g;

				// calculate solution for successive orders
				for (i = 1; i < N; i++)
				{

					e *= (1.0f - g * g);
					if (e == 0.0f)
					{
						throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
					}

					// calculate inner product
					inner = R[i + 1];
					for (j = 0, l = i; j < i; j++, l--)
					{
						inner += a[j] * R[l];
					}

					// update prediction coefficients
					g = -(inner / e);
					for (j = 0, l = i - 1; j < i; j++, l--)
					{
						Z[j] = a[j] + g * a[l];
					}

					// copy vector
					for (j = 0; j < i; j++)
					{
						a[j] = Z[j];
					}

					a[i] = g;

				}
			}

			return a;
		}

		/// <summary>
		/// Invert a symmetric square Toeplitz matrix.
		/// </summary>
		/// <param name="T">The left-most column of the symmetric Toeplitz matrix.</param>
		/// <returns>The inverse matrix.</returns>
		/// <exception cref="ArgumentNullException">
		/// <B>T</B> is a null reference.
		/// </exception>
		/// <exception cref="RankException">
		/// The length of <B>T</B> must be greater than zero.
		/// </exception>
		/// <exception cref="SingularMatrixException">
		/// The Toeplitz matrix or one of the the leading sub-matrices is singular.
		/// </exception>
		/// <remarks>
		/// This static member combines the <b>UDL</b> decomposition and Trench's algorithm into a
		/// single algorithm. When compared to the non-static member it requires minimal data storage
		/// and suffers from no speed penalty.
		/// <para>
		/// Trench's algorithm requires <b>N</b> squared FLOPS, compared to <b>N</b> cubed FLOPS
		/// if we simply solved a linear Toeplitz system with a right-side identity matrix (<b>N</b> is the matrix order).
		/// </para>
		/// </remarks>
		public static FloatMatrix Inverse(FloatVector T)
		{

			FloatMatrix X;

			// check parameters
			if (T == null)
			{
				throw new System.ArgumentNullException("T");
			}
			else if (T.Length < 1)
			{
				throw new System.RankException("The length of T must be greater than zero.");
			}
			else if (T.Length == 1)
			{
				X = new FloatMatrix(1);
				X[0, 0] = 1.0f / T[0];
			}
			else
			{

				int N = T.Length;
				float f, g;
				int i, j, l, k, m, n;
				X = new FloatMatrix(N);

				// calculate the predictor coefficients
				FloatVector Y = FloatSymmetricLevinson.YuleWalker(T);

				// calculate gamma
				f = T[0];
				for (i = 1, j = 0; i < N; i++, j++)
				{
					f += T[i] * Y[j];
				}
				g = 1.0f / f;

				// calculate first row of inverse
				X[0, 0] = g;
				for (i = 1, j = 0; i < N; i++, j++)
				{
					X[0, i] = g * Y[j];
				}

				// calculate successive rows of upper wedge
				for (i = 0, j = 1, k = N - 2; i < N / 2; i++, j++, k--)
				{
					for (l = j, m = i, n = N-1-j; l < N - j; l++, m++, n--)
					{
						X[j, l] = X[i, m] + g * (Y[i] * Y[m] - Y[k] * Y[n]);
					}
				}

				// this is symmetric matrix ...
				for (i = 0; i <= N / 2; i++)
				{
					for (j = i + 1; j < N - i; j++)
					{
						X[j, i] = X[i, j];
					}
				}

				// and a persymmetric matrix.
				for (i = 0, j = N - 1; i < N; i++, j--)
				{
					for (k = 0, l = N - 1; k < j; k++, l--)
					{
						X[l, j] = X[i, k];
					}
				}

			}

			return X;
		}

		#endregion Public Static Methods

	}

}
